package com.sensetime.opencvdemo;

import java.util.logging.Logger;
import no.uib.cipr.matrix.DenseMatrix;
import no.uib.cipr.matrix.EVD;
import no.uib.cipr.matrix.Matrix;
import no.uib.cipr.matrix.NotConvergedException;

/**
 * 
 */
public class PrincipalComponentAnalysis {

    private static Logger logger = Logger
            .getLogger(PrincipalComponentAnalysis.class.getSimpleName());

    private double[][] rawData = null; // 原始数据
    private double[][] principalData = null; // 主成分提取后的数据
    private double[] principalEigenValues = null; // 主要的特征值
    private double[][] principalEigenVectors = null; // 主要的特征向量
    private double[] eigenValues = null; // 特征值
    private double[][] eigenVectors = null; // 特征向量
    private int[] selected = null; // 选择的特征值序号
    private double proportion = 0.9; // 主成分的比重，默认为90%

    /**************** 构造函数里配置一些参数 ****************/
    public PrincipalComponentAnalysis() {

    }

    public PrincipalComponentAnalysis(double proportion) {
        this.proportion = proportion;
    }

    /**************** get方法 ****************/
    public double[][] getRawData() {
        return rawData;
    }

    public double[][] getPrincipalData() {
        return principalData;
    }

    public double[] getPrincipalEigenValues() {
        return principalEigenValues;
    }

    public double[][] getPrincipalEigenVectors() {
        return principalEigenVectors;
    }

    public double[] getEigenValues() {
        return eigenValues;
    }

    public double[][] getEigenVectors() {
        return eigenVectors;
    }

    public int[] getSelected() {
        return selected;
    }

    public double getProportion() {
        return proportion;
    }

    /**************** PCA的内部方法 ****************/
    // 将原始数据标准化
    private double[][] calcStandardlizer(double[][] rawData) {
        double[][] standardData = null;
        if (rawData != null) {
            int N = rawData.length; // 二维矩阵的行数，样本个数
            int p = rawData[0].length; // 二维矩阵的列数，属性个数
            // 每个属性对应的列是该属性的一个采样，近似为该属性的分布
            double[] average = new double[p]; // 每一列的平均值
            double[] var = new double[p]; // 每一列的方差
            standardData = new double[N][p]; // 标准化后的列向量组成的矩阵

            // 取得平均值
            for (int k = 0; k < p; k++) {
                double temp = 0;
                for (int i = 0; i < N; i++) {
                    temp += rawData[i][k];
                }
                average[k] = temp / N;
            }
            // 取得方差
            for (int k = 0; k < p; k++) {
                double temp = 0;
                for (int i = 0; i < N; i++) {
                    temp += Math.pow(rawData[i][k] - average[k], 2);
                }
                var[k] = temp / (N - 1);
            }
            // 取得标准化的矩阵 期望为0，方差为1，简化相关系数的计算公式
            for (int i = 0; i < N; i++) {
                for (int j = 0; j < p; j++) {
                    standardData[i][j] = (double) ((rawData[i][j] - average[j]) / Math
                            .sqrt(var[j]));
                }
            }
        } else {
            logger.info("There is no raw data.");
        }
        return standardData;
    }

    // 计算样本相关系数矩阵 各列之间相互计算（共p列），形成p*p的矩阵
    // 输入为标准化之后的矩阵，利用期望=0，方差=1简化了相关系数的计算
    private double[][] calcCoefficientOfAssociation(double[][] standardData) {
        double[][] assosiationMatrix = null;
        if (standardData != null) {
            int n = standardData.length; // 二维矩阵的行号
            int p = standardData[0].length; // 二维矩阵的列号
            assosiationMatrix = new double[p][p];// 相关系数矩阵
            for (int i = 0; i < p; i++) {
                for (int j = 0; j < p; j++) {
                    double temp = 0;
                    for (int k = 0; k < n; k++) {
                        temp += standardData[k][i] * standardData[k][j];
                    }
                    assosiationMatrix[i][j] = temp / (n - 1);
                }
            }
        }
        return assosiationMatrix;
    }

    // 计算相关系数矩阵的特征值
    private double[] calcEigenValue(double[][] assosiationMatrix) {
        // assosiationMatrix是一个方阵
        double[] eigenValues = null;
        if (assosiationMatrix != null) {
            DenseMatrix Assosiation = new DenseMatrix(assosiationMatrix);
            int len = assosiationMatrix.length;
            EVD evd = new EVD(len);
            try {
                evd.factor(Assosiation);
                eigenValues = evd.getRealEigenvalues();
            } catch (NotConvergedException e) {
                e.printStackTrace();
                logger.info(e.getMessage());
            }
        }
        return eigenValues;
    }

    // 计算相关系数矩阵的特征向量
    private double[][] calcEigenVector(double[][] assosiationMatrix) {
        // assosiationMatrix是一个方阵
        double[][] eigenVectors = null;
        if (assosiationMatrix != null) {
            DenseMatrix Assosiation = new DenseMatrix(assosiationMatrix);
            int len = assosiationMatrix.length;
            EVD evd = new EVD(len);
            try {
                evd.factor(Assosiation);
                DenseMatrix tempMatrix = evd.getLeftEigenvectors();
                eigenVectors = new double[len][len];
                for (int i = 0; i < len; i++) {
                    for (int j = 0; j < len; j++) {
                        eigenVectors[i][j] = tempMatrix.get(i, j);
                    }
                }
            } catch (NotConvergedException e) {
                e.printStackTrace();
            }
        }
        return eigenVectors;
    }

    // 假设阈值是90%，选取最大的前几个特征值的序号
    private int[] selectPrincipalComponent(double[] eigenValues) {
        int p = eigenValues.length; // 特征值个数
        int[] sortedEigenValueIndex = new int[p]; // 特征值由大到小，序号排序
        for (int i = 0; i < p; i++) {
            sortedEigenValueIndex[i] = i;
        }
        double[] temp = new double[p]; // 特征值副本
        System.arraycopy(eigenValues, 0, temp, 0, p);

        // TODO 排序过程可以再优化
        double tmp = 0.0;
        int pos = 0;
        for (int i = 1; i <= p; i++) {
            for (int j = 0; j < p - i - 1; j++) {
                if (temp[j] < temp[j + 1]) {
                    tmp = temp[j];
                    temp[j] = temp[j + 1];
                    temp[j + 1] = tmp;
                    pos = sortedEigenValueIndex[j];
                    sortedEigenValueIndex[j] = sortedEigenValueIndex[j + 1];
                    sortedEigenValueIndex[j + 1] = pos;
                }
            }
        }

        double total = 0.0; // 特征值的和
        for (int i = 0; i < p; i++) {
            total += temp[i];
        }

        int count = 0;
        double sum = 0.0;
        for (int i = 0; i < p; i++) {
            if (sum / total <= proportion) {
                sum += temp[i];
                count++;
            }
        }
        int[] selected = new int[count];
        System.arraycopy(sortedEigenValueIndex, 0, selected, 0, count);
        return selected;
    }

    // 取得主特征值
    private double[] calcPrincipalEigenValues(double[] eigenValues,
                                              int[] selected) {
        int p = eigenValues.length;
        double[] principalEigenValues = new double[selected.length];
        for (int i = 0; i < selected.length; i++) {
            principalEigenValues[i] = eigenValues[selected[i]];
        }
        return principalEigenValues;
    }

    // 取得主特征向量，即变换矩阵
    private double[][] calcPrincipalEigenVectors(double[][] eigenVectors,
                                                 int[] selected) {
        int p = eigenVectors.length;
        double[][] principalEigenVectors = new double[p][selected.length];
        for (int i = 0; i < selected.length; i++) {
            for (int j = 0; j < p; j++) {
                principalEigenVectors[j][i] = eigenVectors[j][selected[i]];
            }
        }
        return principalEigenVectors;
    }

    // 原始数据的主成分数据
    private double[][] calcPrincipalComponent(double[][] rawData) {
        Matrix A = new DenseMatrix(rawData);
        Matrix B = new DenseMatrix(principalEigenVectors);
        Matrix C = new DenseMatrix(rawData.length,
                principalEigenVectors[0].length);
        A.mult(B, C); // C=A*B

        double[][] principalData = new double[C.numRows()][C.numColumns()];
        for (int i = 0; i < C.numRows(); i++) {
            for (int j = 0; j < C.numColumns(); j++) {
                principalData[i][j] = C.get(i, j);
            }
        }
        return principalData;
    }

    /**************** PCA的主流程 ****************/
    public void buildPCA(double[][] rawData) {
        this.rawData = rawData;
        double[][] standardData = calcStandardlizer(rawData);
        double[][] assosiationMatrix = calcCoefficientOfAssociation(standardData);
        this.eigenValues = calcEigenValue(assosiationMatrix);
        this.eigenVectors = calcEigenVector(assosiationMatrix);
        this.selected = selectPrincipalComponent(eigenValues);
        this.principalEigenValues = calcPrincipalEigenValues(eigenValues,
                selected);
        this.principalEigenVectors = calcPrincipalEigenVectors(eigenVectors,
                selected);
        this.principalData = calcPrincipalComponent(rawData);
    }
}
